/***************************************************************************
 *
 *   Copyright (C) 2003 by Haydon Knight
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
#include "dsp/Transformation.h"
#include "dsp/TimeSeries.h"
#include "dsp/BitSeries.h"

#include "dsp/TimeOrder.h"

dsp::TimeOrder::TimeOrder()
 : Transformation<TimeSeries,BitSeries>("TimeOrder",outofplace)
{
  rapid = Polarisation;
  bitseries_ndat = -1;  // ie use the timeseries' ndat
  offset = 0;
}

dsp::TimeOrder::~TimeOrder(){ }

void dsp::TimeOrder::transformation(){
	  
  if(input->get_nbit()!=32)
    throw Error(InvalidState,"dsp::ChannelOrder::transformation()",
		"nbit!=32- this routine is only set up for floats generated by dsp::TimeOrder\n");

  int64_t bs_ndat = 0;

  if( bitseries_ndat<0 )
    bs_ndat = input->get_ndat()-offset;
  else if( bitseries_ndat+offset > input->get_ndat() )
    bs_ndat = input->get_ndat()-offset;
  else
    bs_ndat = bitseries_ndat;

  if( verbose )
    fprintf(stderr,"dsp::TimeOrder::transformation() will output a BitSeries of ndat " I64 " offset by " UI64 " samps (input ndat=" UI64 ")\n",
	    bs_ndat, offset, input->get_ndat());

  if( bs_ndat <= 0 )
    throw Error(InvalidParam,"dsp::TimeOrder::transformation()",
		"Your bitseries ndat came out to be 0 or negative (precisely: " I64 ").  This means your offset (" UI64 ") and your requested bitseries ndat (" UI64 ") were incompatible with the actual timeseries ndat (" UI64 ").  Try reducing your offset.",
		bs_ndat, offset, bitseries_ndat, input->get_ndat());
		
  output->Observation::operator=( *input );
  output->change_start_time( offset );
  output->resize( bs_ndat );

  // number of floats between (t0,f0) and (t1,f0) of a BitSeries
  register const unsigned output_stride = output->get_nchan()*output->get_npol()*output->get_ndim();

  register const unsigned input_stride = input->get_ndim();

  if( verbose )
    fprintf(stderr,"dsp::TimeOrder::transformation() got output_stride=%d input_stride=%d swap=%s\n",
	    output_stride, input_stride, input->get_swap()?"true":"false");

  if( !input->get_swap() ){

    for( unsigned idim=0; idim<input->get_ndim(); idim++){
      for( unsigned ichan=0; ichan<input->get_nchan(); ichan++){
	for( unsigned ipol=0; ipol<input->get_npol(); ipol++){
	  float* in = ((float*)input->get_datptr(ichan,ipol))+idim + input_stride*offset;
	  float* out = 0;
	  if( rapid==Channel )
	    out = (float*)output->get_rawptr() + idim + input->get_ndim()* (ichan+output->get_nchan()*ipol);
	  else if( rapid==Polarisation )
	    out = (float*)output->get_rawptr() + idim + input->get_ndim()* (ipol+output->get_npol()*ichan);

	  if( verbose )
	    fprintf(stderr,"dsp::TimeOrder::transformation(): idim=%d ichan=%d ipol=%d in=%p out = %p+%ld = %p as rapid=%s\n",
		    idim,ichan,ipol,in,output->get_rawptr(),out-(float*)output->get_rawptr(),
		    out, (rapid==Channel)?"Channel":(rapid==Polarisation)?"Polarisation":"unknown");

	  register const unsigned isamp_end = bs_ndat*input->get_ndim(); 
	  register unsigned output_samp=0;
	  
	  if( verbose )
	    fprintf(stderr,"dsp::TimeOrder::transformation() got idim=%d ichan=%d ipol=%d input->data=%p in=%p diff=%ld output->data=%p out=%p diff=%ld\n",
		    idim, ichan, ipol,
		    input->get_datptr(0,0), in, in-input->get_datptr(0,0),
		    output->get_rawptr(), out, out-(float*)output->get_rawptr()); 

	  for( unsigned isamp=0; isamp<isamp_end;
	       isamp+=input_stride, output_samp+=output_stride)
	    out[output_samp] = in[isamp];
	  
	}
      }
    }     
    
    
  }
  
  else {
 
    unsigned idim;

    for( idim=0; idim<input->get_ndim(); idim++){
      for( unsigned ichan=0; ichan<input->get_nchan()/2; ichan++){
	for( unsigned ipol=0; ipol<input->get_npol(); ipol++){
	  float* in = ((float*)input->get_datptr(ichan,ipol))+idim + input_stride*offset;
	  float* out = 0;
	  if( rapid==Channel )
	    out = (float*)output->get_rawptr() + idim + input->get_ndim()* ((ichan+input->get_nchan()/2)+output->get_nchan()*ipol);
	  else if( rapid==Polarisation )
	    out = (float*)output->get_rawptr() + idim + input->get_ndim()* (ipol+output->get_npol()*(ichan+input->get_nchan()/2));
	  
	  register const unsigned isamp_end = bs_ndat*input->get_ndim(); 
	  register unsigned output_samp=0;
	  
	  for( unsigned isamp=0; isamp<isamp_end;
	       isamp+=input_stride, output_samp+=output_stride)
	    out[output_samp] = in[isamp];
	}
      }
    }

    for( idim=0; idim<input->get_ndim(); idim++){
      for( unsigned ichan=input->get_nchan()/2; ichan<input->get_nchan(); ichan++){
	for( unsigned ipol=0; ipol<input->get_npol(); ipol++){
	  float* in = ((float*)input->get_datptr(ichan,ipol))+idim + input_stride*offset;
	  float* out = 0;
	  if( rapid==Channel )
	    out = (float*)output->get_rawptr() + idim + input->get_ndim()* ((ichan-input->get_nchan()/2)+output->get_nchan()*ipol);
	  else if( rapid==Polarisation )
	    out = (float*)output->get_rawptr() + idim + input->get_ndim()* (ipol+output->get_npol()*(ichan-input->get_nchan()/2));	  

	  register const unsigned isamp_end = bs_ndat*input->get_ndim(); 
	  register unsigned output_samp=0;
	  
	  for( unsigned isamp=0; isamp<isamp_end;
	       isamp+=input_stride, output_samp+=output_stride)
	    out[output_samp] = in[isamp];
	}
      }
    }
    
    output->set_swap( false );
  }

  if( verbose )
    fprintf(stderr,"Bye from dsp::TimeOrder::transformation()\n");
}
